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One dimensional S = 1 antiferromagnetic Heisenberg chains with bond alternation and uniaxial 
single-ion-type anisotropy are studied by numerical exact diagonalization of finite size systems. 
We calculate the ground state phase diagram using the twisted boundary condition method for 
D >Q. The ground states of this model contain the Haldane state, large- _D state and the dimer 
state. Using conformal field theory and the level spectroscopy method, we calculate the D and 
(5-dependence of the Luttinger liquid paraments K and the critical exponent v of the energy gap. 
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phase 



§1. Introduction 

One dimensional antiferromagnetic Heisenberg chains 
have been the subject of recent investigations by numer- 

' ous groups. The various phenomena of current interest 
are shown to be due to quantum effects. A uniform anti- 
ferromagnetic Heisenberg chain is known to have a gap- 

\ less ground state for half-integer-spin. Especially, the 

' exact solution is available for S = 1/2 chain. In con- 

\ trast, for the integer-spin,^) there is a gap between the 

■ first excited state and ground state. 

] One dimensional 5 = 1 antiferromagnetic Heisen- 
berg chains with bond alternation undergo a transi- 

\ tion from the Haldane state to the dimer state as the 
strength of bond alternation increases. This transi- 

\ tion has been studied by many methods. On the 

' other hand, the phase boundaries and critical proper- 

] ties of S' = 1 Heisenberg chains with uniaxial single-ion 
anisotropy^' were determined by Glaus and Schneider us- 

] ing a finite size scaling method. The physical picture of 
both states was clarified by Nijs and Rommelse in terms 

\ of string order'^' and by Tasaki using a stochastic geo- 
metric representation.®' Using an exact diagonalization 
method and a phenomenological renormalization group 
approach, Tonegawa et al. calculated the phase diagram 

\ of the S — 1 antiferromagnetic Heisenberg chain with 
bond alternation and uniaxial single-ion-type anisotropy. 
They found the Haldane phase, Neel phase, dimer phase 
and large-D phase in the ground state. No phase tran- 
sition was found, however, between the dimer phase and 
large-D phase. ^' In this paper, we calculate the ground 
state phase diagram of a one demensional S = \ an- 
tiferromagnetic Heisenberg chain with bond alternation 
and uniaxial single-ion-type anisotropy using the twisted 
boundary condition method.^' This improves the ac- 
curacy of the transition points. Using conformal field 
theory and the level spectroscopy method, we obtain the 
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Luttinger liquid parameter K and the critical exponent 
V of the energy gap. 

This paper is organized as follows. In the next section, 
the model Hamiltonian is defined and the numerical re- 
sults are presented. From the exact diagonalization data 
with the twisted boundary condition, the ground state 
phase diagram is obtained. The Luttinger liquid param- 
eter K is calculated by conformal field theory and the 
level spectroscopy method at the transition points. The 
final section is devoted to a summary and discussion. 

§2. Numerical Results 

2.1 Model Hamiltonian 

The Hamiltonian is given by 

N N 

^ = - {-^y5}SiSi+, +Dj2{sn', (2.1) 

1=1 1=1 

where Si is a spin-1 operator. The parameter S{—1 < 
S < 1) and D{D > 0) represent the bond alternation 
and uniaxial single-ion anisotropy, respectively. We take 
(5 > without loss of generality. The ground state of this 
model can be regarded as a dimer state and large-D state 
for large S and D. For D = 0, the ground state is the 
Haldane state or the dimer state depending as to whether 
S < Sc or S > Sc (<5c = 0.260 ± 0.001), respectively. 
For (5 = 0, the ground state is the Haldane state or a 
large-D state according to D < Dc or D > Dq (Dc = 

1.001 ± 0.001), respectively. According to ref. 9, there 
is no transition point between the dimer state and the 
large-D state. 

2. 2 Phase diagram for D > 

In order to determine the phase boundary with high 
accuracy, we use the twisted boundary method of Ki- 
tazawa and Nomura.^' The Hamiltonian is numeri- 
cally diagonalized to calculate the two low lying energy 
levels with the twisted boundary condition {S^_^_^ = 
-Sf, Sl^, = -Si 5^+1 = Sf) for TV = 8,10,12,14 
and 16 by the Lanczos mothed. 



2 



Wei Chen, Kaauo Hida and Bryan C. Sanctuary 



Prom the results of Tonegawa et. al.,^^ it is known 

that two different kinds of ground states are reahzed 
for I? > and for —1 < (5 < 1. For small value of 
6 and D, the ground state is the Haldane state with a 
valence bond solid (VBS) structure. Under the twisted 
boundary condition, the eigenvalues of the space inver- 
sion P {Si S'jv_i+i) and the spin reversal T {Sf 
-S^, -Sf) arc aU equal to -1.5. lo) As 6 and/or 

D increases, a phase transition takes place to the dimer 
state or the large- Z) state for which P = 1 and T — 1. We 
make use of the P and T eigenvalues to distinguish the 
Haldane phase and dimer or large- 1? phase. For example, 
if the value of S is fixed, the energies of the two states vary 
with D. For small D, the energy of the Haldane state 
is lower than that of the dimer or large-D state. As D 
increases, the dimer or large-D state becomes lower than 
the Haldane state. The two levels cross at one point. 
This is the finite size transition point D = Dc{N). Fig- 
ure 1 shows the D-dcpendence of the two lowest levels 
for TV = 16 and S = 0.05. We extrapolate the critical 
point as DdN) = L»c(oo) + aN-"^ + bN-*. Figure 2 
shows the extrapolation procedure. 



-13.6 
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Fig. 1. The D-dcpcndcncc of tlic two lowest energy eigenvalues 
with twist boundary condition. The energies of the Haldane 
state and dimer state or large-D state are represented by o and 
•, respectively, for AT = 16 and 5 = 0.05. 




In Fig. 3, the ground state phase boundary is repre- 
sented by the solid line. Inside the solid line, the ground 
state is the Haldane state. Outside the solid line, the 
ground state is the dimer or the large-D state. We find 
no singularity between the dimer and the large-D state 
and these two states can be regarded as a single phase. 
This is in agreement with ref. 9 

2.3 Critical behaviour on the phase boundary 

At the critical points, we calculate the energy for 
the magnetization = at fixed wave number k by 
the Lanczos exact diagonalization method under peri- 
odic boundary conditions. The magnetization is = 

N 

^Sf. The system sizes are N = 8,10,12,14 and 16. 

1=1 

The ground state has = and k = 0. 

It is known that the finite size correction to the ground 
state energy is related to the central charge c and the spin 
wave velocity vs as follows, ^^"-^^^ 



N 



Eg{N) 



N , 



TTCVs 

6N^' 



^;s= lim —[Ek,iN)-E^ 

N—yoo ZTT 



(2.2) 



(2.3) 



Here Eg{N) is the ground state energy, and Ek^iN) is 
the energy of the excited state with wave number ki = ^ 
and magnetization M^- = 0. Also, c is the central charge 
and Eoo is the ground state energy per unit cell in the 
thermodynamic limit. After appropriate extrapolation 
to ^ oo. we have checked that the central charge c is 
close to unity on the phase boundary as shown in Fig. 4. 
Error bars are estimated from the difference between the 
extrapolated values with the formulae c{N) = c+Ci/N'^ 
and c{N) =c+Ci/N^ + C2/N^. Therefore, we expect 
that the present model can be described by a Gaussian 
model on the critical line, 



Fig. 2. The extrapolation procedure of finite size Dc for (5 = 0.05. 
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Fig. 3. Phase diagram of this model. The Haldane state is inside 
the solid line, the dimer state and large-D state are outside the 
solid line. 



Fig. 4. The D dependence of the numerically obtained central 
charge c. 



dx 



(2.4) 



where (p is the boson field defined in the interval < 
(j) < \/27r, n is the momentum density conjugate to (f> 
which satisfies [(j){x),Il{x')] = iS{x — x') and K is the 
Luttinger liquid parameter. Deviation from the critical 
point is described by the term cos \/20. The low energy 
properties of our model can be described by the one di- 
mensional quantum double sine-Gordon theory near the 
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critical line, 



n 



SG 



+ 



1 

27ra2 



Ax 



(2.5) 



j dx cos \/20 -|- ^^^^ j dx cos \/8^. 



The last term becomes marginal at the SU(2) symmet- 
ric point {D = 0). For _D > 0, it becomes irrelevant. 
We keep this term in the Hamiltonian HgQ however be- 
cause it determines the leading finite size correction to 
the critical dimensions as discussed below. 

The scaling dimensions x„ of the operators are re- 
lated with the energy eigenvalues En{N, M^) of the cor- 
responding excited states = limjv— >c!o Xn{N) where, 



XniN) 



N 
27rfs 



[Er,{N,M')-Eg{N,M')]. (2.6) 



We can identify the correspondence between the op- 
erators in boson representation and the eigenstates of 
spin chains by comparing their symmetry properties.^' 
Let us denote the scaling dimensions of the operators 
•\/2cos-\/2'^ and •\/2sin\/2(?!> by xi and X2, respectively. 
Both should be equal to K/2 in the thermodynamic limit 
as determined from their correlation functions. Numeri- 
cally, these exponents arc calculated using eq. (2.6) for 
= 0, P = 1, A; = 0[a;i] and for = 0, P = -1, 
k = 0[x2]- Thus the value of K can be determined from 
these values. To eliminate the finite size correction due 
to y2, it is more convenient to use the combination, 



K{N)=xi{N)+X2{N), 



(2.7) 



as proposed by Kitazawa and Nomura,^) because the 
leading finite size connections to xi and X2 are given 

by 



Xi 



X2 



2^ 2^N' 



K 

y 



--(-) 

2 ^N' 



2K-2 



(2.8) 



(2.9) 



We assume the formula 
K 



K{N) 



(2.10) 



Ci C2 

and extrapolate K{N) to iV — ^ 00. The extrapolation 
procedure is shown in Fig. 5. The A'^-dependence of 
K{N) is non-monotonic for small 5. Actually it increases 
with N for relatively small N but starts to decrease for 
large N . Although such non-monotonic behavior is not 
observed for large (5, we expect a similar behavior for 
N > 16. Therefore we do not include the data for small 
size systems in the extrapolation procedure of K{N) and 
only use the values for N = 12, 14 and 16. Further, we re- 
gard these extrapolated values as an upper bound while 
the lower bound is estimated from the value for iV = 16 
in the regime 0.15 < i5 < 0.256. The ZJ-dependence of K 
is shown in the Fig. 6. For D>OA, K is larger than 2. 
This implies that the effective interaction between the el- 
ementery fermionic excitation is attractive. In contrast, 
for D<0.4, K is less than 2 and the effective interaction 
is repulsive. 



Since the second term of eq. (2.5) produces the gap, 
the critical exponent v of the energy gap is given by 

' ' (2.11) 



2-xi 



K 




0.005 



0.01 



1/N 



2 0.015 



Fig. 5. The extrapolation procedure for finite size K with D = 
1.001,5 = 0.000; D = 0.985,5 = 0.050; D = 0.931,5 = 
0.100; D = 0.837,5 = 0.150; D = 0.670,5 = 0.200; D = 
0.400, 5 = 0.242; D = 0.200, 5 = 0.256 and D = 0.100, 5 = 0.260 
from top to bottom. 




Fig. 6. The Z?-dependence of K. 



The D dependence of the numerically obtained critical 
exponent u is shown in Fig. 7. For 5 = 0, £' = 1.001 the 
value of is 1.465 ± 0.022. This is consistent with the 
estimation of Glaus and Schneider^) who give u = 1.5, 
Our estimation procedure is more accurate. 

Our model approachs the SU(2) symmetric point in 
the limit D ^ 0. This gives v = 2/3. Our value of v 
is consistent with this and approaches 2/3 in the same 
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Fig. 7. The D-dependence of the numerically obtained critical 
exponent v. 



Fig. 8. The A dependence of the numerically obtained critical 
exponent v of model (3.1) (filled symbols). The open symbols 
are the same exponent for the present model (2.1). 



limit. 

§3. Summary and Discussion 

The ground state phase diagram of a spin-1 antifer- 
romagnotic Heisenberg chain with bond alternation and 
uniaxial single-ion-type anisotropy is accurately deter- 
mined by a numerical diagonalization method which uses 
the twisted boundary condition for I? > 0. The energy 
gap exponent u is calculated by analyzing the numerical 
diagonalization data using conformal field theory. 

The results of this work are compared with our previ- 
ous work^^) on the one dimensional 5=1/2 XXZ model 
with dimerization 1 — j and quadrumerization 5. The 
Hamiltonian is given by 

2N 

^ = '^^j{^2l-1^2l + ^21-1^21 + ^^2l-1^2l) 
1=1 

2N 

+ ^(1 + + SIS^^^, 

1=1 

+ A5|,52Vi), (3.1) 

where A is the anisotropy parameter. Figure 8 shows 

the A and 5 dependences of the critical exponent ly. For 
S < 6c, the critical line of this model tends to that of 
the present model in limit j — > — oo and A — > 1 with 
the product j{A — 1) being kept fixed. The product 
j{A — 1) corresponds to D/2 in the present model. In 
Fig. 8, we plot the exponent v of the model (3.1) against 
A for 5 = 0, 0.2 and 0.3. The open symbols at A = 1 
show the corresponding values of the present model. For 
5 = 0, the exponent v monotonically increases with A 
up to 1.465 ±0.022. For 6 = 0.2 and 0.3, however, a non 
monotonic A-dependence and strong (5-dependence of v 
near the isotropic point is noted. This is due to the fact 
that the critical line of the model (3.1) tends to j — !■ — oo 
in the isotropic limit. For 6 > 6c, it decreases down to 
2/3. 

Our numerical result shows that the value of K varies 



from 1 to 2.633 ± 0.021 as D varies from to 1.001 ± 
0.001. In the neighbourhood oi K = 2 our model can 
be described by a weakly interacting fermion gas. The 
explicit construction of an effective weak coupling theory 
is also a natural extension of our work and is left for 
future studies. 

The numerical calculation is performed using the HI- 
TAC S820 and SR2201 at the Information Processing 
Center of Saitama University and the FACOM VPP500 
at the Supercomputer Center of Institute for Solid State 
Physics, University of Tokyo. This work is supported by 
a Grant-in- Aid for Scientific Research from the Ministry 
of Education, Science, Sports and Culture of Japan and a 
research grant from the Natural Science and Engineering 
Research Council of Canada (NSERC). 
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In the numerical calculation to determine the phase boundary, the D-term was missing 
on the boundary site i = N. Due to the strong system size dependence caused by this error, 
the phase boundary and exponents were inaccurate. Although there are corrections to all 
figures, we only present the corrected Figs. 1, 2, 3, 6 and 7. We have checked that c is closer 
to unity than the data in original Fig. 4. The data points in Fig. 5 move slightly, but we do 
not reproduce it here because Fig. 5 is used only to explain the way of extrapolation. The 
correction to the open symbols in Fig. 8 is unvisibly small. 

In the text, the following numerical values should be corrected: Dc = 0.968 it 0.001 at 
6 = 0. This value of Dc clearly deviates from unity. Therefore the speculation concerning the 
hidden symmetry at D = 1 and 5 = is withdrawn. The exponent ly = 1.4732 it 0.0004 at 
6 = 0. The parameter K varies from 1 to 2.6424 zt 0.0003. The system size in eqs. (2,2), 
(2.3) and (2.6) should read N/2. 

The authors are grateful to Hiroki Nakano and Torn Sakai for pointing out our error. 
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